12 Solid-liquid phase change
In this chapter, we will focus on ‘cold’ phase-change processes, and in particular a continuous material that undergoes melting or solidification. Considering solidifying (freezing) water it is, for instance, of interest to ask for the rate at which the solidification occurs, which translates into the growth rate of the resulting layer of ice with time. This type of problem setting is classically known as the Stefan problem. It is named after the Austrian scientist Josef Stefan (1835-1893), who developed the theory and compared it to data that had been acquired during several polar expeditions (1879-1885) to the North-West Passage. The Stefan problem assumes the phase-change to take place in a homogeneous, continuum material at rest. In reality, however, the situation is often more complicated due to the presence of thermal and compositional convection or mechanical forces. While all of these situations are considerably more complex, the phase-change process as it has originally been described by Stefan is still their fundamental building block. We will hence look at the classical Stefan problem first, and later extend our discussion to complex settings.
12.1 The one-phase Stefan problem
The Stefan problem in 1d aims at determining the temperature evolution \(T(x,t)\) within a in a domain \(x_0 \leq x \leq x_1\) subject to phase-change.
Let’s assume that the domain is initially filled with water ice, hence
\[T(x,0) = T_{ini} \leq T_m,\]
in which \(T_m\) is the melting temperature of ice, hence \(T_{ini}\) is colder than that. We now start heating from the left at location \(x=0\), where we fix the temperature to be \(T_l\). We observe that the ice will start to melt, and that the water-ice interface located at \(X_m(t)\) will propagate through the domain from left to right. Hence, the water portion continuously increases.
The special case, in which the initial ice temperature corresponds to the melting temperature \(T(x,0) = T_m\) is referred to as the isothermal ice case. This is contrasted by a situation, in which the initial temperture is lower than the melting temperature, which is referred to as polythermal ice.
Here, we will consider isothermal initial conditions for. The temperature evolution in the liquid water is governed by the heat equation, which takes the form \[ \begin{align*} \partial_t T = \alpha \, \partial_x^2 T \qquad \text{for} \qquad 0 \leq x < X_m(t). \end{align*} \tag{12.1}\]
As introduced in Chapter 11, \(\alpha\) denotes the thermal diffusivity given by \(\tfrac{\kappa}{\rho c_p}\). The heat equation denotes a parabolic, second-order partial differential equation, for which two boundary conditions are necessary. According to our problem setting, the left domain boundary and the temperature at the water-ice interface are given:
\[ T(x_0,t) = T_l \qquad \text{and} \qquad T(X_m(t),t) = T_m \]
The system is not yet closed, however, as we do not know the position of the interface \(X_m(t)\), which is a function of time. We need an additional piece of information, which is given by the local energy balance at the interface, a heat flux jump condition commonly referred to as the Stefan condition. For the one-phase Stefan problem it reads:
\[ \begin{align*} \underbrace{\rho_s L \underbrace{\partial_t X_m(t)}_{\substack{\text{propagation speed of} \\ \text{the water-ice interface}}}}_{\text{energy needed for the phase-change}} = \underbrace{- \kappa \partial_x T (\underbrace{X_m^{-}(t)}_{\substack{\text{water side of} \\ \text{the interface}}},t)}_{\text{specific heat flux}} \end{align*} \tag{12.2}\]
The Stefan condition includes the latent heat of melting \(L\) \([J/kg]\), which is also referred to as the enthalpy of fusion. For water ice we have \(L \approx 334\) \(J/kg\). This means that the rate at which the interface propagates is proportional to the local heat flux at the interface, hence the heat available to melt the ice. The two equations Equation 12.1 and Equation 12.2 together comprise the so-called ‘one-phase Stefan problem’
To solve the system, we employ the fact that the solution to the heat equation is self-similar. We introduce the similarity variable \(\zeta = x/\sqrt{t}\), and assume the temperature to be a function of this similarity variable only, hence \(T(\zeta)\). The differentials translate according to
\[ \begin{align*} \partial_t T & = \frac{d}{d \zeta} T \frac{d \zeta}{d t} = - \frac{1}{2} \frac{x}{t \sqrt{t}} \frac{d}{d \zeta} T \\ \partial_x T & = \frac{d}{d \zeta} T \frac{d \zeta}{d x} = \frac{1}{\sqrt{t}} \frac{d}{d \zeta} T \\ \partial^2_x T & = \frac{d^2}{d \zeta^2} T \frac{1}{\sqrt{t}} \frac{d \zeta}{d x} = \frac{1}{t} \frac{d^2}{d \zeta^2} T \end{align*} \]
Substitution into the heat Equation 12.1 yields
\[ \begin{align*} \frac{\zeta}{2} \frac{d}{d \zeta} T = - \alpha \frac{d^2}{d \zeta^2} T. \end{align*} \tag{12.3}\]
This is now an ordinary differential equation that depends on \(\zeta\) only, and can easily be solved subject to the specified boundary conditions. Back substitution of \(\zeta = x/\sqrt{t}\) leaves us with
\[ \begin{align*} T(x,t) = T_l - (T_l - T_m) \frac{\text{erf}\left( \frac{x}{2 \sqrt{\alpha t}}\right)}{\text{erf}\left( \frac{X_m}{2 \sqrt{\alpha t}}\right)}, \end{align*} \tag{12.4}\]
in which \(\text{erf}\) denotes the error function
\[\text{erf}(x) := \tfrac{2}{\sqrt{\pi}} \int_0^x e^{-y^2} dy.\]
While we now have an analytic expression for the temperature evolution in the liquid water, we are not yet able to actually evaluate and quantify it, as we do not know the position of the interface \(X_m\), which is a parameter in Equation 12.4.
Considering the previously introduced similarity structure, however, results in the following Ansatz:
\[ \begin{align*} \frac{X_m(t)}{\sqrt{t}} = C, \end{align*} \tag{12.5}\]
in which \(C\) is a constant. Knowing \(C\) would allow us to infer ton the interface position \(X_m\) for a given time. Out of convinience, we factor out \(2 \sqrt{\alpha}\) and write
\[ \begin{align*} X_m(t) = 2 \sqrt{\alpha} \lambda \sqrt{t}, \end{align*} \] such that the unknown coefficent is
\[\lambda = \frac{X_m(t)}{2 \alpha \sqrt{t}}.\]
This doesn’t change the fact that instead of searching for \(X_m\) directly, we reformulated the problem and now look for \(\lambda\), a constant parameter.
Substituting this Ansatz as well as the temperature evolution Equation 12.4 into the Stefan condition Equation 12.2 yields:
\[ \begin{align*} \frac{1}{2} \lambda 2 \sqrt{\alpha} \frac{1}{\sqrt{t}} = \frac{c_p}{L} \alpha \frac{-(T_l-T_m)}{\text{erf}{\tfrac{X_m}{2 \sqrt{\alpha t}}}} \frac{2}{\sqrt{\pi}} e^{-\tfrac{X^2_m}{4 \alpha t}} \frac{1}{2 \sqrt{\alpha t}}, \end{align*} \]
which can be simplified into an algebraic equation for \(\lambda\):
\[ \begin{align*} \lambda \sqrt{\pi} = \underbrace{\frac{c_p (T_m - T_l)}{L}}_{=:Ste} \frac{1}{\text{erf}(\lambda) e^{\lambda^2}} = \frac{Ste}{\text{erf}(\lambda) e^{\lambda^2}} \end{align*} \tag{12.6}\]
Parameter \(Ste\) denotes the the dimensionless Stefan Number, which is the ratio between sensible heat and latent heat. Note, that Equation 12.6 cannot be solved explicitly for \(\lambda\). Instead, the function \(g\) deinfed according to
\[ g(\lambda) := Ste - \underbrace{\sqrt{\pi} \lambda e^{\lambda^2} \text{erf}\lambda}_{=G(\lambda)} = Ste - G(\lambda) = 0 \] is solved for its (unique) root, for instance, via a Newton iteration. Once the root has been determine, we can determine the position of the phase-change interface according to Equation 12.5. After that, the temperature field Equation 12.5 can be evaluated, granting us the full solution.
12.2 The two-phase Stefan problem
A similar solution strategy can be constructed for the two-phase Stefan problem, in which the initial temperature differs from the melting temperature (see exercises).
12.3 Numerical solution techniques
Before discussing more complex physical situations, we will discuss potential numerical solution strategies. In general, we distinguish between interface tracking and fixed-grid methods. Here, we will discuss a fixed-grid method, referred to as the enthalpy method:
In order to derive the enthalpy formulation, we consider a small reference volume \(V\) that consists of a solid part \(V_s\) and a liquid part \(V_l\). The respective volume fractions are defined as
\[\phi_i = \frac{|V_i|}{|V|}, \quad i \in {l,s}.\]
We consider a fully saturated situation, in which \(\phi_l + \phi_s = 1\). We also assume the volume to be in local thermal equilibrium, such that we do not have to distinguish between solid and liquid temperature \(T_l = T_s = T\). The governing component-wise energy balance in enthalpy form can then be derived using a local volume average mixture model, e.g. following Bear (1972):
\[ \begin{align*} \partial_t \left( \phi_i H_i \right) + \nabla \cdot \left(\phi_i H_i \mathbf{v}_i \right) - \nabla \cdot \left( \phi_i \kappa_i \nabla T \right) + \text{interface terms} = 0 \end{align*} \tag{12.7}\]
with \(i \in \{s,l\}\) and interface terms of similar magnitude but opposite sign. The phase enthalpies are defined as, see also Voller, Swaminathan, and Thomas (1990):
\[ \begin{align*} H_s(\mathbf x, t) = \int_{T_{ref}}^{T} \rho_s c_{p,s} d \theta \quad \text{and} \quad H_l(\mathbf x, t) = \int_{T_{ref}}^{T} \rho_l c_{p,l} d \theta + \rho_l L. \end{align*} \]
Note, that the enthalpy is defined with respect to a chosen reference temperature \(T_{ref}\). Adding up solid and liquid enthalpy equations yields an equation for the mixture enthalpy \(H\)
\[ \begin{align*} \partial_t H + \nabla \cdot \left(\phi_s H_s \mathbf{v}_s + \phi_l H_l \mathbf{v}_l \right) = \nabla \cdot \left( \kappa \nabla T \right), \end{align*} \]
in which \(H\) and \(\kappa\) are given by
\[ \begin{align*} H = \phi_s \int_{T_{ref}}^{T} \rho_s c_{p,s} d \theta + \phi_l \int_{T_{ref}}^{T} \rho_l c_{p,l} d \theta + \phi_l \rho_l L \quad \text{and} \quad \kappa = \phi_s \kappa_s + \phi_l \kappa_l. \end{align*} \]
In the absence of convection we get
\[ \begin{align*} \partial_t H = \nabla \cdot \left( \kappa \nabla T \right), \end{align*} \]
which again seems to pose a closure problem, as \(H\) and \(T\) are both unknowns in the system. However, enthalpy and temperature are no independent quantities, we can understand the equation either as a solution in terms of \(T\) or a solution in terms of \(H\). Let’s consider option 1 first. We find
\[ \partial_t H(T) = \underbrace{\phi_s \rho_s c_{p,s} + \phi_l \rho_l c_{p,l} + \frac{d}{dT}\left(\phi_l \rho_l L\right)}_{=:\rho c_{p,a}} \, \partial_t T, \]
in which \(\rho c_{p,a}\) is non-linear in \(T\) and referred to as the apparent heat capacity. The equation to be solved reads
\[ \begin{align*} \rho c_{p,a} (T) \partial_t T = \nabla \cdot \left( \kappa \nabla T \right). \end{align*} \]
The apparent heat capacity is parametrized as a function of temperature. This approach to the numerical solution of phase-change problems is hence referred to as the apparent heat capacity method. Alternatively, one can parametrize the porosity as a function of temperature. We consider this situation for the special case \(\rho = \rho_s = \rho_l\) and \(c_p = c_{ps} = c_{pl}\), which yields
\[ H = \rho c_{p} (T-T_{ref}) + \phi \rho L. \]
Substitution into the energy equation in enthalpy form yields
\[ \partial_t ( \rho c_{p} (T-T_{ref}) + \phi \rho L ) = \rho c_p \partial_t T + \rho L \partial_t \phi = \nabla \cdot \left( \kappa \nabla T \right). \]
Again we have multiple options to proceed. Either, we parametrize \(\phi(T)\) as a heavy side function at melting temperature \(T_m\), which implies that \(\tfrac{d}{dT} \phi\) is given as a Kronecker \(\delta\). This results in an apparent heat capacity method via identifying
\[ \rho c_{p,a} : = \rho c_p (1 + \tfrac{L}{c_p}\delta(T_m)). \]
Alternatively, we can transfer the phase-change related term to the right hand side and parametrize directly. This reads
\[ \rho c_p \partial_t T = \nabla \cdot \left( \kappa \nabla T \right) - \rho L \partial_t \phi, \]
and is commonly referred to as a source term method.
Finally, we can also solve the energy equation directly in \(H\), which in contrast to the previous approaches requires us to parametrize \(H(T)\). From a mathematical well-posedness point of view, this even seems to be the method of choice, as we avoid a parametrization of strong non-linerities, and yields an invertible expression, such that it makes sense to write
\[ \partial_t H = \nabla \cdot \left( \kappa \nabla H^{-1} \right), \]
which is commonly referred to as enthalpy method.
Enthalpy-type methods are contrasted by interface-tracking methods, which are inspired by the semi-analytical Stefan problem. These integrate the heat equation on either side of the phase interface and subsequently evaluate the local heat flux jump (Stefan condition) to infer on the local interface velocity. Interface tracking methods tend to be very precise, yet are technically demanding. They are appropriate only if the physical situation, we are considering corresponds to a distinct interface, meaning that the solid and liquid phase are separated by a clearly defined interface that can be modeled by a jump discontinuity (example: freezing fresh water). However, other cold phase interface morphologies are possible also, e.g.:
- Mushy layer: As a result of solidification, solid and liquid phase are separated by a morphologically complex ‘porous’ layer that arises from dendritic growth of the solid material. The liquid phase might get trapped in liquid lenses within the solid (example: aqueous solutions or metal alloys)
- Continuous interface: Solid and liquid phase are separated by a transition regions in which the phase’s volume fractions gradually change from one pure phase region to the other (example: glass)
Both mushy layers, and continuous transition regions can be perceived as continuous interfaces in phase-change, and are better modeled by the previously introduced fixed-grid methods.